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Abstract - Large scale simulations and analytical theory have been combined to obtain the 
non-equilibrium velocity distribution, f{v), of randomly accelerated particles in suspension. The 
simulations are based on an event-driven algorithm, generalised to include friction. They reveal 
strongly anomalous but largely universal distributions which are independent of volume fraction 
and collision processes, which suggests a one-particle model should capture all the essential fea- 
tures. We have formulated this one-particle model and solved it analytically in the limit of strong 
damping, where we flnd that f{v) decays as 1/v for multiple decades, eventually crossing over to 
a Gaussian decay for the largest velocities. Many particle simulations and numerical solution of 
the one-particle model agree for all values of the damping. 



Introduction — In recent years there has been growing in- 
L^^jt?erest in so called active matter, referring to the ability of 
tlie constituents to move actively by either extracting en- 
ergy from the environment or depleting an internal energy 
\^^46pot. Examples are motor proteins, bacterial swimmers 
^\0T motile cells [T]. Whereas the mechanism that drives 
Q^t|he individual active particle has been studied for many 
years [2H1] , the collective behavior of a large number of in- 
ff^ dividuals has been addressed only recently. Very rich be- 
^^havior has been observed, ranging from pattern formation 
and nonequilibrium phase transitions to turbulence [51|5] . 
' . I Active particles on mesoscopic to macroscopic scales have 
I>also been realized in the form of self-propelled colloids 
(Janus particles) [7] and vibrated polar granular rods [5]. 
^ I More generally, granular particles that are driven by ran- 
Ci dom kicks may be considered active matter with, however, 
the direction of motion being random. 

Our focus here are the velocity distributions of active 
particles in suspension. Whereas in equilibrium, the ve- 
locities universally follow the Maxwell-Boltzmann distri- 
bution, this does not hold for nonequilibrium stationary 
states, where in general deviations from the Maxwell- 
Boltzmann distribution are observed. Few studies have 
focused on the velocity distribution in the context of ac- 
tive cell and bacteria suspensions [Hllin]- In M exten- 
sive experimental data were taken for several cell types. 
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allowing for a statistical analysis of the cell's velocities. 
The authors concluded that exponential distributions are 
a general characteristic feature of cell motility. Such ex- 
ponential distributions have indeed been found in models 
of active Brownian particles jllj : however other distribu- 
tions have been seen as well, depending on the mechanism 
of self-propulsion [Tn[T2] . 

For driven granular media on the other side, numer- 
ous studies have been performed to analyze velocity dis- 
tributions. In experiments, various driving mechanisms 
were shown to produce non-Gaussian velocity distribu- 
tions [TSUlSj . If the particle's motion is strongly damped 
either due to the surrounding fluid or due to collisions 
with the wall, the velocity distributions are exponential. 
In [16| the authors use a single-particle simulation of a 
frictional particle to explain the observed velocity distri- 
bution. Their argument was turned into a Fokker-Planck 
equation [111 [5D], whose stationary solution is in good 
agreement with experiment |19| . 

In this Letter we study a simple model of active particles 
in a suspension, described below, using event-driven simu- 
lations. We obtain nearly universal velocity distributions 
which depend primarily on a single parameter, and which 
exhibit significant deviations from Gaussian behavior, but 
also non-exponential tails (see Fig. [T]). Further, we de- 
velop a single-particle theory that shows good agreement 
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with the simulation data. 

Model — Here we discuss a simple model for active parti- 
cles: hard spheres placed in a fluid with a viscous drag 7, 
that are accelerated at discrete times and undergo elastic 
collisions. 

The equation of motion for particle i reads 



Av, 



At 



Av,: 



coll 



At 



(1) 



Dr 



The driving force is modeled as discrete kicks with am- 
plitude Ap = mAv and frequency for- The components 
of the kick velocity, e.g. Av^, are drawn from a Gaussian 
distribution with mean and variance : 
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and for the other components accordingly. We ignore hy- 
drodynamic interactions. 

The above dynamics is a very crude approximation to 
the run-and-tumble behavior of bacteria, such as E.Coli 
and others pTH2^ . In a time interval At, a particle is ac- 
celerated once and subsequently performs a random mo- 
tion determined by the surrounding fluid and interactions 
with the other particles. If the bacteria acceleration events 
(strokes) are sufficiently rare, subsequent kicks may be re- 
garded as are uncorrelatcd in direction. 

We are interested in a steady state, where the energy 
due to dissipation is balanced by the energy input due to 
random kicks: 



(3) 



where d is the dimensionality of the system. In the follow- 
ing wc will choose units such that lengths are measured 
in units of particle radius and mass in units of particle 
mass. We choose the time scale so that the average steady 
state kinetic energy is d/2, which corresponds to fcsT = I 
for a thermal system. In these units the driving ampli- 
tude becomes — 2j/fur, leaving three independent 
parameters: 7, fon and the volume fraction 77. Wc will 
consider moderately dilute systems for which the particle 
collision frequency is well-described by the Enskog result 
'-^coU = ^'^XV/V^ with the Carnahan-Starling expression 
for the pair correlation at contact X = (1 ^ — v)^- 

Thus our three parameters provide three independent time 
scales: 7, //j,., and uJcoii (in place of ij). 

Simulations — We performed event driven simulations of 
hard spheres. The original algorithm [24j[25] was changed 
in order to implement friction as in |26j . The main effort of 
an event-driven simulation of ballistically moving particles 
goes into the calculation whether two particles will collide 
or not. If a collision between particle i and j will occur, 
the difference of their trajectories, 

r,(t) - r,(t) = r,,,(t) = r,j(to) + v,;j(to)(t - to) (4) 

must be equal to the sum of their radii at time tcou, i.e.. 
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Fig. 1: (colour on-line) Velocity distributions for volume frac- 
tion 77 = 0.35, for = ^cou, and several values of /? = 
0.1,1,3,5,10. The dashed-dotted line shows the Maxwell- 
Boltzmann distribution. The coloured solid lines show the first 
iterative solutions of the one-particle model (see text below) for 
^ = 3,5,10. 



yielding a quadratic equation in tcoii — to- For the damped 
motion, 7 7^ 0, one can still integrate the equations of 
motion in between collisions analytically: 

1 _ g-7(*-to) 

W = r,,j(to) + v,,j(to) (6) 

7 

Compared to ballistic motion, the linear time interval 
between two collisions [tcoii — to) is simply replaced by 
(1 — e"'''(*'^°"~*°))/7. Since we know the collision time from 
the ballistic simulation, we can just use the above relation 
to determine the collision times for the damped system. 
The remaining events in the simulation — driving events, 
wall collisions, sub-box wall collisions — are handled ac- 
cordingly. The only remaining difference in the damped 
system is that the place of a collision with another parti- 
cle or a (sub-box) wall might not be within range of the 
damped motion. If this is the case, the collision will not 
occur, instead the particle will slow down until a driving 
event takes place. 

We have simulated a 3-dimensional system of 2122416 
monodisperse spheres with volume fractions 77 = 0.05 and 
0.35. The system is equilibrated with 7 = and no forc- 
ing. Subsequently, damping and the acceleration force are 
switched on. Then, after another 100 collisions per particle 
to ensure relaxation to a stationary state, the velocity dis- 
tribution is measured. These simulations were conducted 
for various values of the parameters 7, for, and 77, with 
the results described as follows. 

In Fig. [1] we show the velocity distribution for rj = 0.35 
and for ~ ^coU with various values of the damping con- 
stant 7. The curves are labeled by the ratio /3 := ^1 for- 
Whereas for very small /3 the distribution is approximately 
Gaussian, we observe increasingly strong deviations for 



p-2 



Anomalous velocity distributions in active Brownian suspensions 




Fig. 2: (colour on-line) Testing the dependence of the velocity 
distribution on the volume fraction rj. Data for rj = 0.05 and 
Tj — 0.35 are shown for 13 — 3 and for /3 = 10. In both cases 
we find little to no dependence on the volume fraction. The 
driving frequency is taken to be for ~ ^^coii- 

larger /3. Small velocities are highly ovcrpopulated with 
an indication of a singularity in the limit of large /3. High 
velocities are overpopulated as well as compared to the 
equilibrium Maxwell-Boltzmann distribution. These devi- 
ations can be understood intuitively as follows: particles 
which have not been recently kicked are damped to nearly 
zero velocity, whereas the recently kicked particles popu- 
late the tail. 

Next, we demonstrate the universality of these distri- 
butions. The three-dimensional parameter space can be 
spanned by the parameters rj, /?, and for/'-^coii- In Fig. [5] 
we test the volume fraction dependence of the velocity dis- 
tribution. For a given value of /?, e.g., /3 = 3, we simulate 
for rj = 0.05 and -q = 0.35, and find no discernable differ- 
ence between the distributions. This holds for all values 
of /3, as shown for /3 = 3, 5, and 10. 

Then, in Fig. [3] we test the dependence of the velocity 
distribution on the ratio for/'^coU- Data for fDr/i^coii = 
1, 10, and 100 are shown for /3 = 3 and for f3 = 5. For 
a specific value of /3, the curves for different for/'^coii he 
essentially on top of each other. 

We summarize the main results of our simulations: 

• While the model contains three independent time 
scales, 7, for and lOcoii, the distribution is almost 
exclusively determined by the ratio /3 = 7//Dr- 

• As a consequence, the distribution is independent of 
volume fraction (or equivalently collision frequency) 
for the investigated range of rj. 

• The one particle velocity distribution is Gaussian only 
in the limit /3 — >■ 0. The distribution shows increas- 
ingly stronger deviations at small and large veloci- 
ties for increasing (3. 



Fig. 3: (colour on-line) Testing the dependence of the velocity 
distribution on the ratio for/^cou for two values of /3. The 
volume fraction is taken to be 77 = 0.05. 

• The distribution seems to develop a singularity at 
small argument as /? — > 00. 

These observations, in particular the insensitivity to col- 
lision rate, have led us to derive an approximate analyti- 
cal theory for the velocity distribution based on a single- 
particle model that neglects collisions. 

Single-Particle-Model — For simplicity, we consider one 
spatial dimension only, assuming that the cartesian com- 
ponents of the velocity are independent. With At = 
1//Dr, we consider the time interval [0, At), within which 
each particle gets one velocity kick at some random time. 
The idea of the calculation is the following: we use the one 
particle distribution at the beginning of the interval as in- 
put and compute the resulting one particle distribution at 
end of the interval, and then require the two distributions 
to be the same in the stationary state. The speed of a 
single particle decreases in At due to damping and gener- 
ally increases due to a velocity kick, denoted by Aw. The 
kick occurs at time r with probability w{t) = provided 
< r < At. We are interested in the velocity distribu- 
tion at the end of the time interval, when the kick velocity 
has decayed to Au/ = Af exp(— 7(At — r)). For a given 
(fixed) kick size Au, this quantity is a random variable 
due to the stochastic occurrence of the kick in the given 
time interval. The conditional probability to find a veloc- 
ity Awy for a given kick size Aw is easily computed from 
the distribution w{t): 

{^■nr— r, e"^ < Awf/Aw < 1 
/3|Aw/| (7) 
else 

To obtain the non-conditional probability, pfe(Aw/) we 
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write 



dAv pkiAvf\Av) P{Av) 



1 1 



PAvf 



Avf 



dAv ■ P{Ai 



(8) 



where P{Av) is the probabihty distribution for the kick 
velocity, given by Eq. ^ with standard deviation a = 

The total velocity at the end of the time interval, 
Vf = Avf+v is the sum of two terms: the kick velocity and 
the velocity from the start of the interval, Vi^ propagated in 
time to the end of the interval, v = Vie~^ . Given the dis- 
tribution of the initial velocities fi{vi), the distribution of 
final velocities (without kick) is given by f{v) = fi{ve^)e^ . 
Since the two velocity contributions Avf and v are sta- 
tistically independent, the probability distribution of the 
sum is given by the convolution: /(w/) = {pk * f)(,Vf). In 
the stationary state, we require that the initial velocity 
distribution is equal to the final velocity distribution. 



du pk{v - u) f{ue'^) 



(9) 



The probability distribution within this single- particle 
model is a function of /3 = ^ / for only, which matches the 
behavior of the many-particle simulation data. With the 
Fourier transform /(fc) = / dv e^^"" f {v) , the above equa- 
tion simplifies. 



and is solved by 



f{k)^Pk{k)f{ke- 



with 



Pk{k) 



dw exp 
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For a given /3, the infinite product can be truncated for 
some value of j ^ 1/(3. 

We now analyze the behavior of this formal solution, 
Eq. ([TT|). in the limits of large and small /3, where we 
can obtain simple analytic expressions for f{v), and for 
intermediate values of (3, where we obtain the distribution 
through an iterative numerical method. 

First, in the /3 — > limit, the Avf distribution pk{Avf ) 
goes to P{Avj), which is a Gaussian. Thus, according 
to Eq. ([S]), the velocity distribution f(v) must map to it- 
self under a convolution with a Gaussian, which requires 
that f{v) must itself be a Gaussian. This stationary limit 
corresponds to a continuous Ornstein-Uhlenbeck process 
[27] . The cumulant relation [55] implied by Eq. (|TT]) speci- 
fies that the variance of f{v) goes to unity as /3 — ?> 0. This 
is confirmed by the simulations with /3 = 0.1, shown in 

Fig.d] 
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Fig. 4: (colour on-line) The three asymptotic solutions from 
Eq. (|13|l (dashed lines) and simulation data for jS — 10, and 
P — 5 in the inset (both simulation data for rj = 0.35). The 
dotted lines depict the range limits from Eq. (|13p . The dashed 
lines are the analytic results from Eq. (|13|l , without any fitting. 



Second, in the large /3 limit, only the j = term in 
Eq. pT]) contributes to the product, and so f{v) =pk{v). 
As such, from Eq. ([S]), we can identify three regions: 




\v\ ^ ere 



(13) 



e-""/4^ \v\ > a 



The middle case corresponds to taking the integration 
range m Eq. (m to be zero to infinity (for positive Avf). 
The top case corresponds to a smooth cutoff to the l/|Aw/| 
behavior as Avf — > 0. The large Avf limit is obtained by 
setting the upper integration limit to infinity, giving the 
complementary error function. In Fig. |3]we plot the veloc- 
ity distribution data for /3 = 10 and /? = 5 and compare to 
the analytic expressions (dashed lines) from Eq. (fT3]) and 
their ranges (dotted lines). The three regions are clearly 
distinguishable and match the simulation data well. Note 
that the l/\v\ region shrinks as (3 decreases. 

Third, for intermediate values of (3 we solved the defin- 
ing Eq. (|9]) numerically by iteration, starting from a 
Maxwell-Boltzmann distribution. The convergence of the 
iteration process is very fast; there is almost no difference 
visible between the first three iterations (see Fig. [5]). To 
quantify the difference between two subsequent iterations, 
we compute the £^ norm of A/(a;) = /"+^(a;) — /"(x). As 
an example, for /3 = 3 we find values of 0(10^'^) between 
the first and second iteration, and 0(10"^) between the 
second and third iteration, respectively. The iterative so- 
lution of Eq. ^ is compared to the data from simulations 
for several values of /? in Fig. [T] No deviations can be 
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Fig. 5: (colour on-line) Main part: The first three iterations 
for (3 — 3 are almost indistinguishable and agree with the sim- 
ulation data. Inset: First iteration and simulation data for 
P = 10, and Pk{v) which is indistinguishable from the first 
iteration. 



detected within the scatter of the data, 
good agreement for all values of /3. 



We find similar 



Conclusion — We have shown that a Brownian suspen- 
sion of interacting particles, subjected to random accel- 
erations, exhibits strongly anomalous velocity distribu- 
tions. An event driven algorithm was generalised to fi- 
nite friction, allowing for large scale simulations of over 
2 million particles. The simulations reveal velocity dis- 
tributions which are universal in the sense that they are 
largely independent of volume fraction and collisions be- 
tween the particles, and only depend on damping rate and 
kick frequency through the ratio /? = for- This has 
led us to consider a simplified one particle model allow- 
ing for an analytical theory of the velocity distribution, 
f{v). The resulting integral equation reduces trivially to 
the Maxwell-Boltzmann distribution for P ^ Q. For large 
/3, we find a divergent distribution for small argument, 
/(O) ^ , a \/v decay for intermediate v and Gaussian 
behavior for the largest argument. Hence there are no 
exponential tails. In Refs. [TTl[2n] an exponential tail was 
obtained for a damped particle kicked by white shot noise, 
but in these works the kick size distribution was exponen- 
tial, rather than the Gaussian we use. For intermediate /3, 
the integral equation for f[x) is solved by iteration with 
very fast convergence. For all /3 we find excellent agree- 
ment between the one particle theory and the simulations. 

Our approach can be generalised in several ways. Both 
the simulations as well as the analytical theory can be gen- 
eralised to other distributions for the kick amplitudes and 
times. Furthermore we plan to study directed motion, po- 
lar particles and rotational degrees of freedom, modeling 
other swimmers. 
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